library(ggplot2)
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/bayes_pvalue_beta0/")

sel=list.files(pattern = "bayes_p.txt")
group1=c("X1T","X2B","M7","M8","M5","M6","M1","M2","M47","M48","M27","M28","M29","M30","M25","M26","M36","M35","M17","M18","M19","M20","M21","M22","M39","M40","M49","M50")
disease=c(rep("SZ",4),rep("BD",6),"SZ","SZ",rep("BD",6),rep("NC",8),rep("SZ",2))
group2=c(rep(c("DCD","DCC"),5),rep("CC",8),rep("NC",8),rep(c("DCD","DCC"),1))

sample_info=data.frame(group1,disease,group2)

rt=data.frame(matrix(NA,1,3,byrow = T))
names(rt)=c("unitID","beta0","group1")
rt=rt[-1,]

for (i in 1:length(sel)) {
file=read.table(sel[i],head=T,sep="\t")
file$unitID=paste(file$chrom,file$position,sep=":")
file$normal_bayes_fdr=p.adjust(file$normal_bayes_pvalue)
file$tumor_bayes_fdr=p.adjust(file$tumor_bayes_pvalue)
#file$num=as.numeric(file$normal_reads1)+as.numeric(file$normal_reads2)+as.numeric(file$tumor_reads1)+as.numeric(file$tumor_reads2)
#file=file[file$num>=10,]
normal=file[file$normal_bayes_pvalue<0.05,]
tumor=file[file$tumor_bayes_pvalue<0.05,]

normal=data.frame(unitID=normal$unitID,beta0=normal$normal_bayes_beta0,group1=unlist(strsplit(gsub(".bayes_p.txt","",sel[i]),"_"))[1])
tumor=data.frame(unitID=tumor$unitID,beta0=tumor$tumor_bayes_beta0,group1=unlist(strsplit(gsub(".bayes_p.txt","",sel[i]),"_"))[2])
rt=rbind(rt,normal,tumor)
}
rt=merge(rt,sample_info,by="group1")

factor1=factor(rt$group1,levels = c("X1T","M7","M5","M1","M47","M49","X2B","M8","M6","M2","M48","M50","M27","M28","M29","M30","M25","M26","M36","M35","M17","M18","M19","M20","M21","M22","M39","M40"))
factor2=factor(rt$group2,levels = c("DCD","DCC","CC","NC"))

ggplot(rt, aes(x=factor1, y=beta0),color=group1) +
  geom_boxplot(aes(fill=factor(group2)))+theme_bw() +
  theme(axis.text.x=element_text(angle=50,hjust=0.5, vjust=0.5)) +
  theme(legend.position="none")

ggplot(rt, aes(x=factor1, y=beta0),color=group1) +
  geom_violin(aes(fill=factor(group2))) +
  theme(axis.text.x=element_text(angle=50,hjust=0.5, vjust=0.5)) +
  theme(legend.position="none")

num1=c(median(rt[rt$group1=="X1T",]$beta0),median(rt[rt$group1=="M7",]$beta0),median(rt[rt$group1=="M5",]$beta0),median(rt[rt$group1=="M1",]$beta0),median(rt[rt$group1=="M47",]$beta0),median(rt[rt$group1=="M49",]$beta0))
num2=c(median(rt[rt$group1=="X2B",]$beta0),median(rt[rt$group1=="M8",]$beta0),median(rt[rt$group1=="M6",]$beta0),median(rt[rt$group1=="M2",]$beta0),median(rt[rt$group1=="M48",]$beta0),median(rt[rt$group1=="M50",]$beta0))
wilcox.test(num1,num2,exact = FALSE)

num1=c(median(rt[rt$group1=="M27",]$beta0),median(rt[rt$group1=="M28",]$beta0),median(rt[rt$group1=="M29",]$beta0),median(rt[rt$group1=="M30",]$beta0),median(rt[rt$group1=="M25",]$beta0),median(rt[rt$group1=="M26",]$beta0),median(rt[rt$group1=="M36",]$beta0),median(rt[rt$group1=="M35",]$beta0))
num2=c(median(rt[rt$group1=="M17",]$beta0),median(rt[rt$group1=="M18",]$beta0),median(rt[rt$group1=="M19",]$beta0),median(rt[rt$group1=="M20",]$beta0),median(rt[rt$group1=="M21",]$beta0),median(rt[rt$group1=="M22",]$beta0),median(rt[rt$group1=="M39",]$beta0),median(rt[rt$group1=="M40",]$beta0))
wilcox.test(num1,num2,exact = FALSE)

twins=c("DCC","NC")
test=rt[rt$group2 %in% twins,]
wilcox.test(beta0 ~ group2, data = test,exact = FALSE)
